Model annual dissimilarity with temperature, fishing, and survey identity

Script 8 for Kitchel et al. 2024 in prep taxonomic diversity manuscript.

library(data.table)
library(MuMIn)
library(ggplot2)
library(cowplot)
library(lme4)
library(stringr)

###Predicts annual dissimilarity with annual characteristics, temperature and fishing values

Pull in - region areas (if not already loaded) - region characteristics (if not already loaded); saveRDS(FishGlob_richness_year_survey, file = here::here(“output”,“FishGlob_richness_year_survey.Rds”)) - fishing (if not already loaded) - temp (if not already loaded)

#physical area by year
region_area_byyear <- fread(here::here("output","region_area_byyear.csv"))

#merged fishing, temp, dissimilarities
dissimilarities_temp_fishing <- fread(here::here("output","dissimilarities_temp_fishing.csv"))

#combine
dissimilarities_temp_fishing_area <- dissimilarities_temp_fishing[region_area_byyear, on = c("survey_unit","year")]

#only jaccard for these analyses
dissimilarities_temp_fishing_area.jaccard <- dissimilarities_temp_fishing_area[dissimilarity_metric == "jaccard_dissimilarity_total",]

Add in average Julian day

#load up julian days
dates_regions <- readRDS(here::here("output","dates_regions.rds"))

#simplify
dates_regions.annual <- unique(dates_regions[,.(survey_unit, year, avg_julian_annual)])

#merge
dissimilarities_temp_fishing_area.jaccard <- dates_regions.annual[dissimilarities_temp_fishing_area.jaccard, on = c("survey_unit","year")]

Pull in palette and name helper

source(here::here("analysis_code","color_links.R"))

Pull in observed trend values

jaccard_total_coefs.r <- fread(here::here("output","jaccard_total_coefs.r.csv"))

Plot fishing and temperature vs. time for all regions

#######TEMPERATURE

#set order by survey unit for plotting
all_surveys <- levels(as.factor(dissimilarities_temp_fishing_area.jaccard$survey_unit))
setorder(dissimilarities_temp_fishing_area.jaccard, survey_unit)

dissimilarities_temp_fishing_area.jaccard[,Survey_Name_Season:=factor(Survey_Name_Season, levels = unique(dissimilarities_temp_fishing_area.jaccard$Survey_Name_Season), ordered = T)]

(sbt_time_survey_facet_1_20 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[1:20] & year > 1979]) +
  labs(y = "Mean bottom temperature (˚C)",  x = "Year") +
  geom_point(aes(y = as.numeric(yearly_mean_bypoint_avg), x = year), alpha = 0.3) +
  geom_smooth(aes(y = as.numeric(yearly_mean_bypoint_avg), x = year), method = "lm") +
  scale_x_continuous(breaks = ~ axisTicks(., log = FALSE)) +
  theme_classic() +
  theme(axis.text.x = element_text(size = 7)) +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4))

ggsave(sbt_time_survey_facet_1_20, path = here::here("figures"), filename = "temp_time_survey_facet_1_20.jpg", height = 12, width =9)

(sbt_time_survey_facet_21_34 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[21:34] & year > 1979]) +
  labs(y = "Mean bottom temperature (˚C)",  x = "Year") +
  geom_point(aes(y = as.numeric(yearly_mean_bypoint_avg), x = year), alpha = 0.3) +
  geom_smooth(aes(y = as.numeric(yearly_mean_bypoint_avg), x = year), method = "lm") +
  scale_x_continuous(breaks = ~ axisTicks(., log = FALSE)) +
  theme_classic() +
  theme(axis.text.x = element_text(size = 7)) +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4))

ggsave(sbt_time_survey_facet_21_34, path = here::here("figures"), filename = "sbt_time_survey_facet_21_34.jpg", height = 12, width =9)

#######FISHING

dissimilarities_temp_fishing_area.jaccard.cc <- dissimilarities_temp_fishing_area.jaccard[complete.cases(dissimilarities_temp_fishing_area.jaccard[,summed_tonnes_scaled_byreg]),]

(fishing_time_survey_facet_1_20 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard.cc[survey_unit %in% all_surveys[1:20] & year > 1979]) +
  labs(y = "Relative fishing catch",  x = "Year") +
  geom_point(aes(y = summed_tonnes_scaled_byreg, x = year), alpha = 0.3) +
  geom_smooth(aes(y = summed_tonnes_scaled_byreg, x = year), method = "lm") +
  scale_x_continuous(breaks = ~ axisTicks(., log = FALSE)) +
  theme_classic() +
  theme(axis.text.x = element_text(size = 7)) +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4))

ggsave(fishing_time_survey_facet_1_20, path = here::here("figures"), filename = "fishing_time_survey_facet_1_20.jpg", height = 12, width =9)

(fishing_time_survey_facet_21_34 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard.cc[survey_unit %in% all_surveys[c(21:34)] & year > 1979]) +
  labs(y = "Relative fishing catch",  x = "Year") +
  geom_point(aes(y = summed_tonnes_scaled_byreg, x = year), alpha = 0.3) +
  geom_smooth(aes(y = summed_tonnes_scaled_byreg, x = year), method = "lm") +
  scale_x_continuous(breaks = ~ axisTicks(., log = FALSE)) +
  theme_classic() +
  theme(axis.text.x = element_text(size = 7)) +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4))

ggsave(fishing_time_survey_facet_21_34, path = here::here("figures"), filename = "fishing_time_survey_facet_21_34.jpg", height = 12, width =9)

NA
NA

Plot fishing and temperature vs. dissimilarity for all regions

#####MEAN TEMP
(preds_sbt_mean_temp_survey_facet_1_20 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[1:20] & year > 1979]) +
  labs(x = "Mean bottom temperature (˚C)",  y = "β-diversity") +
  geom_point(aes(x = as.numeric(yearly_mean_bypoint_avg), y = annual_dissimilarity_value), alpha = 0.3) +
  geom_smooth(aes(x = as.numeric(yearly_mean_bypoint_avg), y = annual_dissimilarity_value), method = "lm") +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
  theme_classic())

ggsave(preds_sbt_mean_temp_survey_facet_1_20, path = here::here("figures"), filename = "preds_sbt_mean_temp_survey_facet_1_20.jpg", height = 12, width =9)

(preds_sbt_mean_temp_survey_facet_21_34 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[21:34] & year > 1979]) +
  labs(x = "Mean bottom temperature (˚C)",  y = "β-diversity") +
  geom_point(aes(x = as.numeric(yearly_mean_bypoint_avg), y = annual_dissimilarity_value), alpha = 0.3) +
  geom_smooth(aes(x = as.numeric(yearly_mean_bypoint_avg), y = annual_dissimilarity_value), method = "lm") +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
  theme_classic())

ggsave(preds_sbt_mean_temp_survey_facet_21_34, path = here::here("figures"), filename = "preds_sbt_mean_temp_survey_facet_21_34.jpg", height = 12, width =9)

#####MINIMUM TEMP
(preds_sbt_min_temp_survey_facet_1_20 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[1:20] & year > 1979]) +
  labs(x = "Minimum bottom temperature (˚C)",  y = "β-diversity") +
  geom_point(aes(x = as.numeric(yearly_min_bypoint_avg), y = annual_dissimilarity_value), alpha = 0.3) +
  geom_smooth(aes(x = as.numeric(yearly_min_bypoint_avg), y = annual_dissimilarity_value), method = "lm") +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
  theme_classic())

ggsave(preds_sbt_min_temp_survey_facet_1_20, path = here::here("figures"), filename = "preds_sbt_min_temp_survey_facet_1_20.jpg", height = 12, width =9)

(preds_sbt_min_temp_survey_facet_21_34 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[21:34] & year > 1979]) +
  labs(x = "Minimum bottom temperature (˚C)",  y = "β-diversity") +
  geom_point(aes(x = as.numeric(yearly_min_bypoint_avg), y = annual_dissimilarity_value), alpha = 0.3) +
  geom_smooth(aes(x = as.numeric(yearly_min_bypoint_avg), y = annual_dissimilarity_value), method = "lm") +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
  theme_classic())

ggsave(preds_sbt_min_temp_survey_facet_21_34, path = here::here("figures"), filename = "preds_sbt_min_temp_survey_facet_21_34.jpg", height = 12, width =9)

#######FISHING

dissimilarities_temp_fishing_area.jaccard.cc <- dissimilarities_temp_fishing_area.jaccard[complete.cases(dissimilarities_temp_fishing_area.jaccard[,summed_tonnes_scaled_byreg]),]

(preds_sbt_mean_fishing_survey_facet_1_20 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard.cc[survey_unit %in% all_surveys[1:20] & year > 1979]) +
  labs(x = "Relative fishing catch",  y = "β-diversity") +
  geom_point(aes(x = summed_tonnes_scaled_byreg, y = annual_dissimilarity_value), alpha = 0.3) +
  geom_smooth(aes(x = summed_tonnes_scaled_byreg, y = annual_dissimilarity_value), method = "lm") +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
  theme_classic())

ggsave(preds_sbt_mean_fishing_survey_facet_1_20, path = here::here("figures"), filename = "preds_sbt_mean_fishing_survey_facet_1_20.jpg", height = 12, width =9)

(preds_sbt_mean_fishing_survey_facet_21_34 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard.cc[survey_unit %in% all_surveys[c(21:34)] & year > 1979]) +
  labs(x = "Relative fishing catch",  y = "β-diversity") +
  geom_point(aes(x = summed_tonnes_scaled_byreg, y = annual_dissimilarity_value), alpha = 0.3) +
  geom_smooth(aes(x = summed_tonnes_scaled_byreg, y = annual_dissimilarity_value), method = "lm") +
  facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
  theme_classic())

ggsave(preds_sbt_mean_fishing_survey_facet_21_34, path = here::here("figures"), filename = "preds_sbt_mean_fishing_survey_facet_21_34.jpg", height = 12, width =9)

NA
NA

###Plot number of tows per year per region

tows_year <- unique(dissimilarities_temp_fishing_area.jaccard[,.(survey_unit, haul_id_count_annual, area_km, year, Survey_Name_Season)])

min_haul_density <- min(tows_year$haul_id_count_annual/tows_year$area_km)
max_haul_density <- max(tows_year$haul_id_count_annual/tows_year$area_km)

(tow_density_survey_facet_1_20 <- ggplot(data = tows_year[survey_unit %in% all_surveys[1:20]]) +
  labs(x = "Year",  y = expression("Tow density (tows per km"^2*")")) +
  geom_point(aes(x = year, y = haul_id_count_annual/area_km)) +
  ylim(c(min_haul_density-0.0001, max_haul_density+0.0001)) +
  facet_wrap(~Survey_Name_Season, ncol = 4) +
  theme_classic())

ggsave(tow_density_survey_facet_1_20, path = here::here("figures"), filename = "tow_density_survey_facet_1_20.jpg", height = 12, width =9)

(tow_density_survey_facet_21_34 <- ggplot(data = tows_year[survey_unit %in% all_surveys[21:34]]) +
  labs(x = "Year",  y = expression("Tow density (tows per km"^2*")")) +
  geom_point(aes(x = year, y = haul_id_count_annual/area_km)) +
  ylim(c(min_haul_density-0.0001, max_haul_density+0.0001)) +
  facet_wrap(~Survey_Name_Season, ncol = 4) +
  theme_classic())

ggsave(tow_density_survey_facet_21_34, path = here::here("figures"), filename = "tow_density_survey_facet_21_34.jpg", height = 12, width =9)

#minimum and maximum dissimilarity value
min_dissimilarity <- min(dissimilarities_temp_fishing_area.jaccard$annual_dissimilarity_value)
max_dissimilarity <- max(dissimilarities_temp_fishing_area.jaccard$annual_dissimilarity_value)

#how does tow density vary with dissimilarity?
tow_density_dissimilarity_1_20 <- 
  ggplot(data = dissimilarities_temp_fishing_area.jaccard[dissimilarity_metric == "jaccard_dissimilarity_total" & survey_unit %in% all_surveys[1:20]]) +
  geom_point(aes(x = (haul_id_count_annual/area_km), y = annual_dissimilarity_value, color = year), size = 2) +
  viridis::scale_color_viridis() +
  facet_wrap(~Survey_Name_Season, ncol = 4) +
  labs(x = "Tow density (tows per km^2)", y = "β-diversity", color = "Year") +
    xlim(c(min_haul_density-0.0001, max_haul_density+0.0001)) +
    ylim(c(min_dissimilarity-0.01,max_dissimilarity+0.01)) +
  theme_classic()

ggsave(tow_density_dissimilarity_1_20, path = here::here("figures"), filename = "tow_density_dissimilarity_1_20.jpg", height = 12, width =9)

tow_density_dissimilarity_21_34 <- 
  ggplot(data = dissimilarities_temp_fishing_area.jaccard[dissimilarity_metric == "jaccard_dissimilarity_total" & survey_unit %in% all_surveys[21:34]]) +
  geom_point(aes(x = (haul_id_count_annual/area_km), y = annual_dissimilarity_value, color = year), size = 2) +
  viridis::scale_color_viridis() +
  facet_wrap(~Survey_Name_Season, ncol = 4) +
  labs(x = "Tow density (tows per km^2)", y = "β-diversity", color = "Year") +
    xlim(c(min_haul_density-0.0001, max_haul_density+0.0001)) +
    ylim(c(min_dissimilarity-0.01,max_dissimilarity+0.01)) +
  theme_classic()

ggsave(tow_density_dissimilarity_21_34, path = here::here("figures"), filename = "tow_density_dissimilarity_21_34.jpg", height = 12, width =9)

###Set up dredge to identify best performing models

options(na.action = "na.fail")

dissimilarity_covariates_dredge.dt <- dissimilarities_temp_fishing_area.jaccard[,.
                  (year, survey_unit,
                    yearly_mean_bypoint_avg, yearly_max_bypoint_avg, yearly_min_bypoint_avg,yearly_seas_bypoint_avg,
                    yearly_mean_bypoint_SD, yearly_max_bypoint_SD, yearly_min_bypoint_SD,yearly_seas_bypoint_SD,
                    yearly_mean_bypoint_avg.s, yearly_max_bypoint_avg.s, yearly_min_bypoint_avg.s,yearly_seas_bypoint_avg.s,
                    annual_dissimilarity_value,
                    haul_id_count_annual,
                    spp_count_annual, depth_annual_avg,
                    depth_annual_range, latitude_annual_avg,
                    latitude_annual_range, area_km, avg_julian_annual, summed_tonnes_scaled_byreg)]

#merge in with colors for plotting predictions by survey
dissimilarity_covariates_dredge.dt <- color_link[dissimilarity_covariates_dredge.dt, on = "survey_unit"]

#If NA for any covariate, delete row
View(dissimilarity_covariates_dredge.dt)
#Deleted:
  #Before 1980 and after 2019
  #Gulf of Saint Laurence South (no depth data)
  #No clear SAU match for Rockall Plateau
dissimilarity_covariates_dredge.dt <- dissimilarity_covariates_dredge.dt[complete.cases(dissimilarity_covariates_dredge.dt)]


dissimilarity_covariates_dredge.dt[, yearly_mean_bypoint_avg.scaledacrossall := scale(yearly_mean_bypoint_avg)][, yearly_mean_bypoint_SD.scaledacrossall := scale(yearly_mean_bypoint_SD)][, yearly_min_bypoint_avg.scaledacrossall := scale(yearly_min_bypoint_avg)][, yearly_max_bypoint_avg.scaledacrossall := scale(yearly_max_bypoint_avg)][, yearly_seas_bypoint_avg.scaledacrossall := scale(yearly_seas_bypoint_avg)][,haul_id_count_annual.scaledacrossall := scale(haul_id_count_annual)][,spp_count_annual.scaledacrossall := scale(spp_count_annual)][,depth_annual_avg.scaledacrossall := scale(depth_annual_avg)][,depth_annual_range.scaledacrossall := scale(depth_annual_range)] [,latitude_annual_avg.scaledacrossall := scale(latitude_annual_avg)][,latitude_annual_range.scaledacrossall := scale(latitude_annual_range)][,area_km.scaledacrossall := scale(area_km)][,julian_scaled := scale(avg_julian_annual)]

###Full model

temp + survey_unit + fishing + area + latitude range + latitude average + depth range + depth average + spp count + # of hauls + julian day + AR for year

For temperature, we will look at: -mean (scaled across all regions) -max (scaled across all regions) -min (scaled across all regions) -seas (scaled across all regions) -SD

Comparing temp variables AND the presence/absence of a temporal autocorrelation term

global_mod_mean_temp_gls <- gls(annual_dissimilarity_value ~ 
                             survey_unit*yearly_mean_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                             julian_scaled,
                           correlation = corAR1(form=~year | survey_unit),
                             data = dissimilarity_covariates_dredge.dt)

global_mod_mean_temp_lm <- lm(annual_dissimilarity_value ~ 
                             survey_unit*yearly_mean_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                             julian_scaled,
                             data = dissimilarity_covariates_dredge.dt)

global_mod_max_temp_lm <- lm(annual_dissimilarity_value ~ 
                             survey_unit*yearly_max_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                            julian_scaled,
                             data = dissimilarity_covariates_dredge.dt)

global_mod_max_temp_gls <-gls(annual_dissimilarity_value ~ 
                             survey_unit*yearly_max_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                            julian_scaled,
                            correlation = corAR1(form=~year | survey_unit),
                             data = dissimilarity_covariates_dredge.dt)

global_mod_min_temp_lm <- lm(annual_dissimilarity_value ~ 
                             survey_unit*yearly_min_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                            julian_scaled,
                             data = dissimilarity_covariates_dredge.dt)

global_mod_min_temp_gls <-gls(annual_dissimilarity_value ~ 
                             survey_unit*yearly_min_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                            julian_scaled,
                            correlation = corAR1(form=~year | survey_unit),
                             data = dissimilarity_covariates_dredge.dt)

global_mod_seas_temp_lm <- lm(annual_dissimilarity_value ~ 
                             survey_unit*yearly_seas_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                             julian_scaled,
                             data = dissimilarity_covariates_dredge.dt)

global_mod_seas_temp_gls <- gls(annual_dissimilarity_value ~ 
                             survey_unit*yearly_seas_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                             julian_scaled,
                             correlation = corAR1(form=~year | survey_unit),
                             data = dissimilarity_covariates_dredge.dt)

global_mod_SD_temp_lm <- lm(annual_dissimilarity_value ~ 
                             survey_unit*yearly_mean_bypoint_SD.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                           julian_scaled,
                             data = dissimilarity_covariates_dredge.dt)

global_mod_SD_temp_gls <- gls(annual_dissimilarity_value ~ 
                             survey_unit*yearly_mean_bypoint_SD.scaledacrossall + #temp and survey unit (possible interaction)
                             survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
                             area_km.scaledacrossall + #area
                             latitude_annual_range.scaledacrossall + #latitude range
                             latitude_annual_avg.scaledacrossall + #latitude avg
                             depth_annual_range.scaledacrossall + #depth range
                             depth_annual_avg.scaledacrossall + #depth avg
                             spp_count_annual.scaledacrossall + #spp #
                             haul_id_count_annual.scaledacrossall +
                           julian_scaled,
                              correlation = corAR1(form=~year | survey_unit),
                             data = dissimilarity_covariates_dredge.dt)

View(AICc(global_mod_mean_temp_lm, global_mod_max_temp_lm, global_mod_min_temp_lm, global_mod_seas_temp_lm, global_mod_SD_temp_lm,
          global_mod_mean_temp_gls, global_mod_max_temp_gls, global_mod_min_temp_gls, global_mod_seas_temp_gls, global_mod_SD_temp_gls))

#build data table to report AICc
global_mod_temp_table <- data.table(
  `Temporal autocorrelation` = c(
    c(rep(F,5),rep(T,5))
  ),
  `Temperature variable` = rep(c(                  
              "Average mean SBT",
              "Average maximum SBT",                        
              "Average minimum SBT",                  
              "Average SBT seasonality",                 
              "SBT SD"),2),
  deltaAICc = signif(min(AICc(global_mod_mean_temp_lm, global_mod_max_temp_lm, global_mod_min_temp_lm, global_mod_seas_temp_lm, global_mod_SD_temp_lm,
          global_mod_mean_temp_gls, global_mod_max_temp_gls, global_mod_min_temp_gls, global_mod_seas_temp_gls, global_mod_SD_temp_gls)[,2])-AICc(global_mod_mean_temp_lm, global_mod_max_temp_lm, global_mod_min_temp_lm, global_mod_seas_temp_lm, global_mod_SD_temp_lm,
          global_mod_mean_temp_gls, global_mod_max_temp_gls, global_mod_min_temp_gls, global_mod_seas_temp_gls, global_mod_SD_temp_gls)[,2],2))

#order by aicc
setorder(global_mod_temp_table,cols = -"deltaAICc")

global_mod_temp_table[,Rank := seq(1,10,by = 1)]

global_mod_sbt_table <- global_mod_temp_table[,.(Rank,`Temporal autocorrelation`,`Temperature variable`, deltaAICc)]

fwrite(global_mod_sbt_table, here::here("output","global_mod_sbt_table.csv"))

Best performing global model includes minimum temperature (centered and scaled) (Still the case as of May 6, 2024)

Now, look at different combinations of all predictors (min temp) using dredge

Global model: global_mod_min_temp

options(na.action = "na.fail") #  prevent fitting sub-models to different datasets
dd <- dredge(global_mod_min_temp)

Fixed term is "(Intercept)"

dd.dt <- data.table(dd)
View(dd)

#only models less than 4 delta AICc (2 models)
dd.dt.2 <- dd.dt[delta <= 2,]

colnames(dd.dt.2)

 [1] "(Intercept)"                                        "area_km.scaledacrossall"                           
 [3] "depth_annual_avg.scaledacrossall"                   "depth_annual_range.scaledacrossall"                
 [5] "haul_id_count_annual.scaledacrossall"               "julian_scaled"                                     
 [7] "latitude_annual_avg.scaledacrossall"                "latitude_annual_range.scaledacrossall"             
 [9] "spp_count_annual.scaledacrossall"                   "summed_tonnes_scaled_byreg"                        
[11] "survey_unit"                                        "yearly_min_bypoint_avg.scaledacrossall"            
[13] "summed_tonnes_scaled_byreg:survey_unit"             "survey_unit:yearly_min_bypoint_avg.scaledacrossall"
[15] "df"                                                 "logLik"                                            
[17] "AICc"                                               "delta"                                             
[19] "weight"                                            

#in this step, we delete coefficient values because we will pull them back in later when we calculate both SE and coefficients (other than interaction, which we keep here)
dd.dt.2.formatted <- dd.dt.2[,Rank := as.numeric(rownames(dd.dt.2))][,.(Rank,
                                `summed_tonnes_scaled_byreg:survey_unit`,
                                `survey_unit:yearly_min_bypoint_avg.scaledacrossall`,
                                delta,
                                weight, survey_unit
                                )]

#add r-squared values
# Iterate through the best models and extract R-squared values
r_squared_values <-list()
for (i in 1:nrow(dd.dt.2.formatted)) {
  summary_data <- summary(get.models(dd, subset = i)[[1]])
  r_squared <- signif(summary_data$r.squared,2)
  r_squared_values <- unlist(c(r_squared_values,r_squared))
}

dd.dt.2.formatted[,"R squared" := r_squared_values]

#empty data table

model_coef_se_fill <- data.table(Rank = as.numeric(), coef_name = as.character(),coef = as.numeric(), se = as.numeric())

for (i in 1:nrow(dd.dt.2.formatted)){
  model_coef_se_single <- data.table(unlist(coefTable(dd,full = T)[[i]]))
  model_coef_se_single[,coef_names := rownames(unlist(coefTable(dd,full = T)[[i]]))]
  model_coef_se_single[,Rank := i]
  
  colnames(model_coef_se_single) <- c("coef","se","df","coef_name","Rank")
  
  #reduce to columns we need
  model_coef_se_single <- model_coef_se_single[,.(Rank, coef_name, coef, se)]
  
  model_coef_se_fill <- rbind(model_coef_se_fill,model_coef_se_single)
}

#format to merge with model rankings and averaged model
model_coef_se_fill[,coef_se := paste0(round(coef,3)," ± ",round(se,3))]

#delete extra columns
model_coef_se_fill <- model_coef_se_fill[,.(Rank,coef_name,coef_se)]

#long to wide
model_coef_se_fill.w <- dcast(model_coef_se_fill, formula = Rank ~ coef_name, value.var = c("coef_se"))

#merge model_coef_se_fill with dd.dt.2.formatted
model_coef_se_AIC <- dd.dt.2.formatted[model_coef_se_fill.w, on = "Rank"]

#model average all models with delta < 4 (50 models)
model_avg_delta4 <-model.avg(dd, subset = delta < 4, fit = T) #NB: The ‘subset’ (or ‘conditional’) average only   averages over the models where the parameter appears. An alternative, the ‘full’ average assumes that a variable is included in every model, but in some models the corresponding coefficient (and its respective variance) is set to zero. Unlike the ‘subset average’, it does not have a tendency of biasing the value away from zero. The ‘full’ average is a type of shrinkage estimator, and for variables with a weak relationship to the response it is smaller than ‘subset’ estimators., fit = T fits the component models again

model_avg_values <- as.data.table(coefTable(model_avg_delta4,fill = T)) # with SE
coef_names <- names(coef(model.avg(dd, subset = delta < 4)))
model_avg_values[,coef_name:=coef_names][,coef:=Estimate][,Estimate:=NULL][,df:=NULL][,se:= `Std. Error`][,`Std. Error` := NULL]

#new column with coef and SE
model_avg_values[,coef_se := paste0(round(coef,3)," ± ",round(se,3))]

#long to wide for model avg
model_avg.wide <- dcast(model_avg_values, formula = . ~ coef_name, value.var = c("coef_se"))

#add rank of "model avg" to table
model_avg.wide[,Rank := "Model avg"]

best_model_sbt_jaccard_table_formatted <- rbind(model_coef_se_AIC, model_avg.wide, fill = T)

#get rid of interaction coefficients
best_model_sbt_jaccard_table_formatted.r <- best_model_sbt_jaccard_table_formatted[,.(Rank,`(Intercept)`,area_km.scaledacrossall,
                                depth_annual_avg.scaledacrossall,
                                depth_annual_range.scaledacrossall,
                                haul_id_count_annual.scaledacrossall,
                                julian_scaled,
                                latitude_annual_avg.scaledacrossall,
                                latitude_annual_range.scaledacrossall,
                                spp_count_annual.scaledacrossall,
                                summed_tonnes_scaled_byreg,
                                survey_unit,
                                yearly_min_bypoint_avg.scaledacrossall,
                                `summed_tonnes_scaled_byreg:survey_unit`,
                                `survey_unit:yearly_min_bypoint_avg.scaledacrossall`,
                                `R squared`,
                                delta,
                                weight
                                )]

#round to 2 significant figures

#names of numeric columns
numeric_cols <- names(best_model_sbt_jaccard_table_formatted.r)[sapply(best_model_sbt_jaccard_table_formatted.r, is.numeric)]

# Apply signif() only to numeric columns
best_model_sbt_jaccard_table_formatted.r[, (numeric_cols) := lapply(.SD, function(x) if (is.numeric(x)) signif(x, digits = 2) else x), .SDcols = numeric_cols]

#change column names, in caption note that all variables are centered and scaled
colnames(best_model_sbt_jaccard_table_formatted.r) <- c(
"Rank",                                              
"Intercept",                                     
"Area",       #scaled across all                    
"Average depth",                  
"Depth range",                
"Number of tows",     
"Average Julian day",
"Average latitude",               
"Latitude range",             
"Species count",                  
"Relative catch", #scaled within region                       
"Survey",                                       
"Average minimum temperature",            
"Survey * relative catch",            
"Survey * avg min temperature",
"R squared",
paste0("\u0394"," AICc"),                                             
paste0("\u03C9"))



#save as csv

fwrite(best_model_sbt_jaccard_table_formatted.r,here::here("output","best_model_sbt_jaccard_table_formatted.csv"))

Take aways Dissimilarity increases with increasing avg # tows Dissimilarity decreases with increasing avg latitude Dissimilarity decreases with increasing Julian day of sampling Dissimilarity increases with increasing latitude range Dissimilarity increases with increasing fishing effort Dissimilarity decreases with increasing depth range Dissimilarity increases with increasing average depth Dissimilarity increases with increasing area Dissimilarity decreases with increasing minimum temperature Dissimilarity increases with increasing # spp Dissimilarity increases with relative fishing catch

Predict dissimilarity across years using averaged model (model_avg_delta4)

dissimilarity_covariates_dredge.dt_predictions <- copy(dissimilarity_covariates_dredge.dt)

  #allowing temp and fishing to vary (aka no changes)
  dissimilarity_covariates_dredge.dt_predictions[,pred_dissim := predict(model_avg_delta4, se.fit = T, full = F)[[1]]][,pred_se := predict(model_avg_delta4, se.fit = T, full = F)[[2]]] #full allows us to switch back to mixed effect models

Warning: argument 'full' ignoredWarning: argument 'full' ignored

  
  
  
  #constant temp in regions (aka take mean of temperature values within survey units so they are the same value within each year for a survey)
  dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg <- copy(dissimilarity_covariates_dredge.dt)
  dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg[,yearly_min_bypoint_avg.scaledacrossall:=mean(yearly_min_bypoint_avg.scaledacrossall),.(survey_unit)]

  #constant fishing in regions (aka take mean of fishing values within survey units so they are the same value within each year for a survey)
  dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg <- copy(dissimilarity_covariates_dredge.dt)
  dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg[,summed_tonnes_scaled_byreg:=mean(summed_tonnes_scaled_byreg),.(survey_unit)]

#and then with consistent temp and fishing in regions (aka take mean of both min temp and fishing)
dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg <- copy(dissimilarity_covariates_dredge.dt)
dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg[,summed_tonnes_scaled_byreg:=mean(summed_tonnes_scaled_byreg),.(survey_unit)][,yearly_min_bypoint_avg.scaledacrossall:=mean(yearly_min_bypoint_avg.scaledacrossall),.(survey_unit)]


#allowing temp and fishing to vary within regs (normal)
    dissimilarity_covariates_dredge.dt_predictions[,pred_dissim := 
                                                     predict(model_avg_delta4, se.fit = T, full = F, newdata =   dissimilarity_covariates_dredge.dt_predictions)[[1]]][,pred_se := 
                                                                                                                                                                         predict(model_avg_delta4, se.fit = T, full = F, newdata =   dissimilarity_covariates_dredge.dt_predictions)[[2]]]

Warning: argument 'full' ignoredWarning: argument 'full' ignored

  
#allowing only fishing to vary within regions (with mean temp)
  
   dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg[,pred_dissim := 
                                                                        predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg)[[1]]][,pred_se :=
                                                                                                                                                                                                              predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg)[[2]]]

Warning: argument 'full' ignoredWarning: argument 'full' ignored

  
  
#allowing only temperature to vary within regions (with mean fishing pressure)
     
   dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg[,pred_dissim :=
                                                                           predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg)[[1]]][,pred_se := 
                                                                                                                                                                                                                                  predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg)[[2]]]

Warning: argument 'full' ignoredWarning: argument 'full' ignored

   
#and then with consistent temp and fishing across regions (mean of both fishing and temp)
dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg <- copy(dissimilarity_covariates_dredge.dt)
dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg[,summed_tonnes_scaled_byreg:=mean(summed_tonnes_scaled_byreg),.(survey_unit)][,yearly_min_bypoint_avg.scaledacrossall:=mean(yearly_min_bypoint_avg.scaledacrossall),.(survey_unit)]

dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg[,pred_dissim := 
                                                                         predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg)[[1]]][,pred_se :=
                                                                                                                                                                                                                 predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg)[[2]]]

Warning: argument 'full' ignoredWarning: argument 'full' ignored

#and then with consistent temp and fishing across regions (mean of both fishing and temp)
dissimilarity_covariates_dredge.dt_predictions_consistenttempfishing <- copy(dissimilarity_covariates_dredge.dt)
dissimilarity_covariates_dredge.dt_predictions_consistenttempfishing[,summed_tonnes_scaled_byreg:=mean(summed_tonnes_scaled_byreg)][,yearly_min_bypoint_avg.scaledacrossall:=mean(yearly_min_bypoint_avg.scaledacrossall)]

dissimilarity_covariates_dredge.dt_predictions_consistenttempfishing[,pred_dissim := 
                                                                         predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistenttempfishing)[[1]]][,pred_se :=
                                                                                                                                                                                                                 predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistenttempfishing)[[2]]]

Warning: argument 'full' ignoredWarning: argument 'full' ignored

#FOR COLOR TO MATCH
#sort color link by survey name season
#alphabetical order
color_link_alpha <- setorder(color_link, survey_unit)

#exclude surveys we don't include
color_link_alpha <- color_link_alpha[survey_unit %in% unique(dissimilarity_covariates_dredge.dt_predictions$survey_unit),]

color_alpha_order <- color_link_alpha[,hex]
label_alpha_order <- color_link_alpha[,Survey_Name_Season]
  

#maintain temp and fishing
#plot
predicted_values_temp_fishing <- ggplot(dissimilarity_covariates_dredge.dt_predictions) +
  geom_point(aes(x = year, y = annual_dissimilarity_value, color = survey_unit)) +
  geom_line(aes(x = year, y = pred_dissim, color = survey_unit)) +
  geom_ribbon(aes(x = year, ymin = pred_dissim-pred_se, ymax = pred_dissim+pred_se, fill = survey_unit), alpha = 0.3) +
  scale_color_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
  scale_fill_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
  labs(x = "Year",y = "Average annual total\nBray Curtis dissimilarity") +
  ylim(0,1.5) +
  theme_classic() +
  ggtitle("Average model predictions")

#average temp and fishing for each region
predicted_values_temp_fishing_meantempfishinginsurvey <- ggplot(dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg) +
  geom_point(aes(x = year, y = annual_dissimilarity_value, color = survey_unit)) +
  geom_line(aes(x = year, y = pred_dissim, color = survey_unit)) +
  geom_ribbon(aes(x = year, ymin = pred_dissim-pred_se, ymax = pred_dissim+pred_se, fill = survey_unit), alpha = 0.3) +
  scale_color_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
  scale_fill_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
  labs(x = "Year",y = "Average annual total\nBray Curtis dissimilarity") +
  ylim(0,1.5) +
  theme_classic() +
  ggtitle("Average model predictions with mean\nsurvey temperature and fishing pressure")

#average temp and fishing across all regions
predicted_values_temp_fishing_meantempfishing <- ggplot(dissimilarity_covariates_dredge.dt_predictions_consistenttempfishing) +
  geom_point(aes(x = year, y = annual_dissimilarity_value, color = survey_unit)) +
  geom_line(aes(x = year, y = pred_dissim, color = survey_unit)) +
  geom_ribbon(aes(x = year, ymin = pred_dissim-pred_se, ymax = pred_dissim+pred_se, fill = survey_unit), alpha = 0.1) +
  scale_color_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
  scale_fill_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
  labs(x = "Year",y = "Average annual total\nBray Curtis dissimilarity") +
  ylim(0,1.5) +
  theme_classic() +
  ggtitle("Average model predictions with mean\novereall temperature and fishing pressure")

#merge plots
predicted_values_sbt_jaccard_fishing_merge <- cowplot::plot_grid(predicted_values_temp_fishing,
                                                 predicted_values_temp_fishing_meantempfishinginsurvey,
                                                 predicted_values_temp_fishing_meantempfishing,
                                                 ncol = 1)

Warning: Removed 27 rows containing missing values (`geom_line()`).Warning: no non-missing arguments to max; returning -InfWarning: no non-missing arguments to max; returning -Inf

ggsave(predicted_values_sbt_jaccard_fishing_merge, path = here::here("figures"), filename = "predicted_values_sbt_jaccard_fishing_merge.jpg", height = 30, width = 14)

Take dissimilarity values from random normal distribution for each year for each region, and then calculate slope (1000 times). Do this for: - Fishing and temperature vary interannually within surveys - Temperature is held constant (as mean over time series for a survey), but fishing varies (allows us to look at relative variance explained) - Fishing is held constant (as mean over time series for a survey), but temperature varies (allows us to look at relative variance explained) - Both fishing and temperature held constant (allows us to see role of other components of the model)

################################################################################
#full predictions
################################################################################
#table with predicted dissimilarity values and standard error of all predicted dissimilarity values (by year)
table <- dissimilarity_covariates_dredge.dt_predictions[,.(survey_unit, pred_dissim, pred_se, year)]
#0) make datatable to populate
  predicted_dissim_trends_rnormruns <- data.table()
#1) NEW PREDICTED VALUES FROM DISTRIBUTION
for (i in 1:1000){
  table[,rnorm_pred := rnorm(1, mean = pred_dissim, sd = pred_se),.(year, survey_unit)]
#2) CALCULATE LINEAR MODEL TO EXTRACT SURVEY SPECIFIC TREND VALUES
  jaccard_total_predicted_lm_singlerun <- lm(rnorm_pred ~ year*survey_unit,data = table)

  model_coefs_reduced_predictions_singlerun <- data.table(summary(jaccard_total_predicted_lm_singlerun)$coefficients)
  model_coefs_reduced_predictions_singlerun[,var := rownames(summary(jaccard_total_predicted_lm_singlerun)$coefficients)]
  
  #limit to interactions only (check this if there are any model changes!) row 2 and rows 34:64
  model_coefs_reduced_predictions_singlerun <- model_coefs_reduced_predictions_singlerun[c(2,34:64),]
  
  #adjust survey unit name by deleting beginning of string
  model_coefs_reduced_predictions_singlerun[,survey_unit := substr(var, 17, str_length(var))][var == "year",survey_unit := "AI"]
  
  #calculate interaction coefficients
  AI_estimate <- model_coefs_reduced_predictions_singlerun[1,Estimate]
  model_coefs_reduced_predictions_singlerun[1,survey_unit_coefficient := AI_estimate]
  model_coefs_reduced_predictions_singlerun[2:32,survey_unit_coefficient := (AI_estimate + Estimate)]
  
  predicted_dissim_trends_rnormruns <- rbind(predicted_dissim_trends_rnormruns, model_coefs_reduced_predictions_singlerun[,.(survey_unit, survey_unit_coefficient)])
  
  print(i)
}
  
#reduce to mean and standard deviation
predicted_dissim_trends_rnormruns[,mean_dissim_coef:= mean(survey_unit_coefficient),survey_unit][,sd_dissim := sd(survey_unit_coefficient),.(survey_unit)]

predicted_dissim_trends_rnormruns.summary <- unique(predicted_dissim_trends_rnormruns[,.(survey_unit, mean_dissim_coef, sd_dissim)])

predicted_dissim_trends_rnormruns.summary[,pred_type := "full"]

################################################################################
#predictions with temp held constant and fishing still varying from year to year
################################################################################
#table with predicted dissimilarity values and standard error of all predicted dissimilarity values (by year)
table_constanttemp <- dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg[,.(survey_unit, pred_dissim, pred_se, year)]
#0) make datatable to populate
  predicted_dissim_trends_rnormruns_constanttemp <- data.table()
#1) NEW PREDICTED VALUES FROM DISTRIBUTION
for (i in 1:1000){
  table_constanttemp[,rnorm_pred := rnorm(1, mean = pred_dissim, sd = pred_se),.(year, survey_unit)]
#2) CALCULATE LINEAR MODEL FOR SLOPE VALUES
  jaccard_total_predicted_lm_singlerun_constanttemp <- lm(rnorm_pred ~ year*survey_unit,data = table_constanttemp)

  model_coefs_reduced_predictions_singlerun_constanttemp <- data.table(summary(jaccard_total_predicted_lm_singlerun_constanttemp)$coefficients)
  model_coefs_reduced_predictions_singlerun_constanttemp[,var := rownames(summary(jaccard_total_predicted_lm_singlerun_constanttemp)$coefficients)]
  
  #limit to interactions only (check this if there are any model changes!) row 2 and rows 34:64
  model_coefs_reduced_predictions_singlerun_constanttemp <- model_coefs_reduced_predictions_singlerun_constanttemp[c(2,34:64),]
  
  #adjust survey unit name by deleting beginning of string
  model_coefs_reduced_predictions_singlerun_constanttemp[,survey_unit := substr(var, 17, str_length(var))][var == "year",survey_unit := "AI"]
  
  #calculate interaction coefficients
  AI_estimate <- model_coefs_reduced_predictions_singlerun_constanttemp[1,Estimate]
  model_coefs_reduced_predictions_singlerun_constanttemp[1,survey_unit_coefficient := AI_estimate]
  model_coefs_reduced_predictions_singlerun_constanttemp[2:32,survey_unit_coefficient := (AI_estimate + Estimate)]
  
  predicted_dissim_trends_rnormruns_constanttemp <- rbind(predicted_dissim_trends_rnormruns_constanttemp, model_coefs_reduced_predictions_singlerun_constanttemp[,.(survey_unit, survey_unit_coefficient)])
  
  print(i)
}
  
#reduce to mean and standard deviation
predicted_dissim_trends_rnormruns_constanttemp[,mean_dissim_coef:= mean(survey_unit_coefficient),survey_unit][,sd_dissim := sd(survey_unit_coefficient),.(survey_unit)]

predicted_dissim_trends_rnormrunsconstant_temp.summary <- unique(predicted_dissim_trends_rnormruns_constanttemp[,.(survey_unit, mean_dissim_coef, sd_dissim)])

predicted_dissim_trends_rnormrunsconstant_temp.summary[,pred_type := "temp_constant"]

predicted_dissim_trends_rnormruns.summary <- rbind(predicted_dissim_trends_rnormruns.summary, predicted_dissim_trends_rnormrunsconstant_temp.summary)

################################################################################
#predictions with fishing held constant (and temperature varying)
################################################################################
#table with predicted dissimilarity values and standard error of all predicted dissimilarity values (by year)
table_constantfishing <- dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg[,.(survey_unit, pred_dissim, pred_se, year)]
#0) make datatable to populate
  predicted_dissim_trends_rnormruns_constantfishing <- data.table()
#1) NEW PREDICTED VALUES FROM DISTRIBUTION
for (i in 1:1000){
  table_constantfishing[,rnorm_pred := rnorm(1, mean = pred_dissim, sd = pred_se),.(year, survey_unit)]
#2) CALCULATE LINEAR MODEL FOR SLOPE VALUES
  jaccard_total_predicted_lm_singlerun_constantfishing <- lm(rnorm_pred ~ year*survey_unit,data = table_constantfishing)

  model_coefs_reduced_predictions_singlerun_constantfishing <- data.table(summary(jaccard_total_predicted_lm_singlerun_constantfishing)$coefficients)
  model_coefs_reduced_predictions_singlerun_constantfishing[,var := rownames(summary(jaccard_total_predicted_lm_singlerun_constantfishing)$coefficients)]
  
  #limit to interactions only (check this if there are any model changes!) row 2 and rows 34:64
  model_coefs_reduced_predictions_singlerun_constantfishing <- model_coefs_reduced_predictions_singlerun_constantfishing[c(2,34:64),]
  
  #adjust survey unit name by deleting beginning of string
  model_coefs_reduced_predictions_singlerun_constantfishing[,survey_unit := substr(var, 17, str_length(var))][var == "year",survey_unit := "AI"]
  
  #calculate interaction coefficients
  AI_estimate <- model_coefs_reduced_predictions_singlerun_constantfishing[1,Estimate]
  model_coefs_reduced_predictions_singlerun_constantfishing[1,survey_unit_coefficient := AI_estimate]
  model_coefs_reduced_predictions_singlerun_constantfishing[2:32,survey_unit_coefficient := (AI_estimate + Estimate)]
  
  predicted_dissim_trends_rnormruns_constantfishing <- rbind(predicted_dissim_trends_rnormruns_constantfishing, model_coefs_reduced_predictions_singlerun_constantfishing[,.(survey_unit, survey_unit_coefficient)])
  
  print(i)
}
  
#reduce to mean and standard deviation
predicted_dissim_trends_rnormruns_constantfishing[,mean_dissim_coef:= mean(survey_unit_coefficient),survey_unit][,sd_dissim := sd(survey_unit_coefficient),.(survey_unit)]

predicted_dissim_trends_rnormrunsconstant_fishing.summary <- unique(predicted_dissim_trends_rnormruns_constantfishing[,.(survey_unit, mean_dissim_coef, sd_dissim)])

predicted_dissim_trends_rnormrunsconstant_fishing.summary[,pred_type := "fishing_constant"]

predicted_dissim_trends_rnormruns.summary <- rbind(predicted_dissim_trends_rnormruns.summary, predicted_dissim_trends_rnormrunsconstant_fishing.summary)


################################################################################
#predictions with both fishing and temperature held constant (variability goes to other factors we don't account for)
################################################################################
#table with predicted dissimilarity values and standard error of all predicted dissimilarity values (by year)
table_constanttempfishing <- dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg[,.(survey_unit, pred_dissim, pred_se, year)]
#0) make datatable to populate
  predicted_dissim_trends_rnormruns_constanttempfishing <- data.table()
#1) NEW PREDICTED VALUES FROM DISTRIBUTION
for (i in 1:1000){
  table_constanttempfishing[,rnorm_pred := rnorm(1, mean = pred_dissim, sd = pred_se),.(year, survey_unit)]
#2) CALCULATE LINEAR MODEL FOR SLOPE VALUES
  jaccard_total_predicted_lm_singlerun_constanttempfishing <- lm(rnorm_pred ~ year*survey_unit,data = table_constanttempfishing)

  model_coefs_reduced_predictions_singlerun_constanttempfishing <- data.table(summary(jaccard_total_predicted_lm_singlerun_constanttempfishing)$coefficients)
  model_coefs_reduced_predictions_singlerun_constanttempfishing[,var := rownames(summary(jaccard_total_predicted_lm_singlerun_constanttempfishing)$coefficients)]
  
  #limit to interactions only (check this if there are any model changes!) row 2 and rows 34:64
  model_coefs_reduced_predictions_singlerun_constanttempfishing <- model_coefs_reduced_predictions_singlerun_constanttempfishing[c(2,34:64),]
  
  #adjust survey unit name by deleting beginning of string
  model_coefs_reduced_predictions_singlerun_constanttempfishing[,survey_unit := substr(var, 17, str_length(var))][var == "year",survey_unit := "AI"]
  
  #calculate interaction coefficients
  AI_estimate <- model_coefs_reduced_predictions_singlerun_constanttempfishing[1,Estimate]
  model_coefs_reduced_predictions_singlerun_constanttempfishing[1,survey_unit_coefficient := AI_estimate]
  model_coefs_reduced_predictions_singlerun_constanttempfishing[2:32,survey_unit_coefficient := (AI_estimate + Estimate)]
  
  predicted_dissim_trends_rnormruns_constanttempfishing <- rbind(predicted_dissim_trends_rnormruns_constanttempfishing, model_coefs_reduced_predictions_singlerun_constanttempfishing[,.(survey_unit, survey_unit_coefficient)])
  
  print(i)
}
  
#reduce to mean and standard deviation
predicted_dissim_trends_rnormruns_constanttempfishing[,mean_dissim_coef:= mean(survey_unit_coefficient),survey_unit][,sd_dissim := sd(survey_unit_coefficient),.(survey_unit)]

predicted_dissim_trends_rnormrunsconstant_tempfishing.summary <- unique(predicted_dissim_trends_rnormruns_constanttempfishing[,.(survey_unit, mean_dissim_coef, sd_dissim)])

predicted_dissim_trends_rnormrunsconstant_tempfishing.summary[,pred_type := "fishing_and_temp_constant"]

predicted_dissim_trends_rnormruns.summary <- rbind(predicted_dissim_trends_rnormruns.summary, predicted_dissim_trends_rnormrunsconstant_tempfishing.summary)

Plotting observed vs predicted

jaccard_fishing_temp_model_observed_predicted_dt <- jaccard_total_coefs.r[predicted_dissim_trends_rnormruns.summary, on = "survey_unit"]

jaccard_fishing_temp_model_observed_predicted_dt[,pred_lower := mean_dissim_coef-sd_dissim][,pred_upper := mean_dissim_coef+sd_dissim]

#FULL MODEL, both temperature and fishing are allowed to vary
jaccard_fishing_temp_model_observed_predicted_lm <- lm(survey_unit_coefficient ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "full"])
summary(jaccard_fishing_temp_model_observed_predicted_lm) #R^2 0.42

Call:
lm(formula = survey_unit_coefficient ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == 
    "full"])

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0040158 -0.0004766 -0.0001640  0.0003970  0.0030020 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -2.510e-06  2.463e-04  -0.010    0.992    
mean_dissim_coef  1.118e+00  2.494e-01   4.483 9.99e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.001385 on 30 degrees of freedom
Multiple R-squared:  0.4011,    Adjusted R-squared:  0.3812 
F-statistic: 20.09 on 1 and 30 DF,  p-value: 9.994e-05

(jaccard_fishing_temp_model_observed_predicted <- ggplot(jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "full"]) +
  geom_errorbar(aes(x = mean_dissim_coef, ymin = lwr, ymax = upr), color = "lightgrey", linewidth = 0.4) +
  geom_errorbarh(aes(y = survey_unit_coefficient, xmin = mean_dissim_coef-sd_dissim, xmax = mean_dissim_coef+sd_dissim), color = "lightgrey", linewidth = 0.4) +
  geom_point(aes(y = survey_unit_coefficient, x = mean_dissim_coef)) +
  geom_smooth(aes(y = survey_unit_coefficient, x = mean_dissim_coef), color = "darkgrey",linetype = "dotted", method = "lm") +
  geom_abline(aes(slope = 1, intercept = 0)) +
  lims(x = c(min(jaccard_fishing_temp_model_observed_predicted_dt$pred_lower),max(jaccard_fishing_temp_model_observed_predicted_dt$pred_upper))) +
  labs(y = "Observed β-diversity trend",x = "Predicted β-diversity trend\n") +
  theme_classic()
)

#fishing constant (fishing constant; temperature varies only)
jaccard_fishing_constant_model_observed_predicted_lm <- lm(survey_unit_coefficient ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "fishing_constant"])
summary(jaccard_fishing_constant_model_observed_predicted_lm) 

Call:
lm(formula = survey_unit_coefficient ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == 
    "fishing_constant"])

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0037145 -0.0009269  0.0001445  0.0012234  0.0026883 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)   
(Intercept)      -0.0002055  0.0002974  -0.691  0.49495   
mean_dissim_coef  1.5387839  0.5107700   3.013  0.00522 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.001568 on 30 degrees of freedom
Multiple R-squared:  0.2323,    Adjusted R-squared:  0.2067 
F-statistic: 9.076 on 1 and 30 DF,  p-value: 0.005221

#Temperature as a predictor, not fishing = R^2 = 0.28 (drop in 12% of variance explained when you lose fishing as predictor)

(jaccard_fishing_constant_model_observed_predicted <- ggplot(jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "fishing_constant"]) +
geom_errorbar(aes(x = mean_dissim_coef, ymin = lwr, ymax = upr), color = "lightgrey", linewidth = 0.4) +
  geom_errorbarh(aes(y = survey_unit_coefficient, xmin = mean_dissim_coef-sd_dissim, xmax = mean_dissim_coef+sd_dissim), color = "lightgrey", linewidth = 0.4) +
  geom_point(aes(y = survey_unit_coefficient, x = mean_dissim_coef)) +
    geom_smooth(aes(y = survey_unit_coefficient, x = mean_dissim_coef), color = "darkgrey",linetype = "dotted", method = "lm") +
  geom_abline(aes(slope = 1, intercept = 0)) +
      lims(x = c(min(jaccard_fishing_temp_model_observed_predicted_dt$pred_lower),max(jaccard_fishing_temp_model_observed_predicted_dt$pred_upper))) +
  labs(y = "Observed β-diversity trend",x = "Predicted β-diversity trend\n(temperature varies fishing constant)") +
  theme_classic()
)

#temp constant (fishing only; temperature constant)
jaccard_temperature_constant_model_observed_predicted_lm <- lm(survey_unit_coefficient ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "temp_constant"])
summary(jaccard_temperature_constant_model_observed_predicted_lm) #0.18 Drop in 22% of variance explained when you lose temperature as a predictor

Call:
lm(formula = survey_unit_coefficient ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == 
    "temp_constant"])

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0037086 -0.0007951 -0.0000619  0.0004231  0.0043027 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)  
(Intercept)      0.0001087  0.0002877   0.378   0.7083  
mean_dissim_coef 0.6701402  0.2673702   2.506   0.0178 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.001627 on 30 degrees of freedom
Multiple R-squared:  0.1731,    Adjusted R-squared:  0.1456 
F-statistic: 6.282 on 1 and 30 DF,  p-value: 0.01785

(jaccard_temperature_constant_model_observed_predicted <- ggplot(jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "temp_constant"]) +
  geom_errorbar(aes(x = mean_dissim_coef, ymin = lwr, ymax = upr), color = "lightgrey", linewidth = 0.4) +
  geom_errorbarh(aes(y = survey_unit_coefficient, xmin = mean_dissim_coef-sd_dissim, xmax = mean_dissim_coef+sd_dissim), color = "lightgrey", linewidth = 0.4) +
  geom_point(aes(y = survey_unit_coefficient, x = mean_dissim_coef)) +
    geom_smooth(aes(y = survey_unit_coefficient, x = mean_dissim_coef), color = "darkgrey",linetype = "dotted", method = "lm") +
  geom_abline(aes(slope = 1, intercept = 0)) +
      lims(x = c(min(jaccard_fishing_temp_model_observed_predicted_dt$pred_lower),max(jaccard_fishing_temp_model_observed_predicted_dt$pred_upper))) +
  labs(y = "Observed β-diversity trend",x = "Predicted β-diversity trend\n(fishing varies temperature constant)") +
  theme_classic()
)

#both temperature and fish held constant
jaccard_fishing_temp_model_observed_predicted_tempfishconstantinsurvey_lm <- lm(survey_unit_coefficient ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "fishing_and_temp_constant"])
summary(jaccard_fishing_temp_model_observed_predicted_tempfishconstantinsurvey_lm) #%11 

Call:
lm(formula = survey_unit_coefficient ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == 
    "fishing_and_temp_constant"])

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0035992 -0.0007889 -0.0002449  0.0009627  0.0039403 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)
(Intercept)      -1.051e-05  3.211e-04  -0.033    0.974
mean_dissim_coef  1.166e+00  8.503e-01   1.372    0.180

Residual standard error: 0.001736 on 30 degrees of freedom
Multiple R-squared:  0.05902,   Adjusted R-squared:  0.02766 
F-statistic: 1.882 on 1 and 30 DF,  p-value: 0.1803

(jaccard_fishing_temp_model_observed_predicted_tempfishconstantinsurvey <- ggplot(jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "fishing_and_temp_constant"]) +
geom_errorbar(aes(x = mean_dissim_coef, ymin = lwr, ymax = upr), color = "lightgrey", linewidth = 0.4) +
  geom_errorbarh(aes(y = survey_unit_coefficient, xmin = mean_dissim_coef-sd_dissim, xmax = mean_dissim_coef+sd_dissim), color = "lightgrey", linewidth = 0.4) +
  geom_point(aes(y = survey_unit_coefficient, x = mean_dissim_coef)) +
    geom_smooth(aes(y = survey_unit_coefficient, x = mean_dissim_coef), color = "darkgrey",linetype = "dotted", method = "lm")+
  geom_abline(aes(slope = 1, intercept = 0)) +
      lims(x = c(min(jaccard_fishing_temp_model_observed_predicted_dt$pred_lower),max(jaccard_fishing_temp_model_observed_predicted_dt$pred_upper))) +
  labs(y = "Observed β-diversity trend",x = "Predicted β-diversity trend\n(fishing and temperature constant)") +
  theme_classic()
)


#merge
jaccard_fishing_sbt_model_observed_predicted_merge <- plot_grid(jaccard_fishing_temp_model_observed_predicted + theme(plot.margin = unit(c(0.1,0.3,0.1,0.1),"cm")),
                                                                     jaccard_fishing_constant_model_observed_predicted + theme(plot.margin = unit(c(0.1,0.3,0.1,0.1),"cm")),
                                                                     jaccard_temperature_constant_model_observed_predicted + theme(plot.margin = unit(c(0.1,0.3,0.1,0.1),"cm")),
                                                                     jaccard_fishing_temp_model_observed_predicted_tempfishconstantinsurvey + theme(plot.margin = unit(c(0.1,0.3,0.1,0.1),"cm")), ncol = 2, labels = c("a.","b.","c.","d."))

`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'

ggsave(jaccard_fishing_sbt_model_observed_predicted_merge, path = here::here("figures"),filename = "jaccard_fishing_sbt_model_observed_predicted_merge.jpg", height =6, width = 8)

NA
NA

Let’s visualize model coefficients with temperature and fishing (similar to figure 2) model_avg_values

model_avg_values

interaction_avg_model_coef <- model_avg_values[c(5,37:99),.(coef_name, coef, se)]

#insert AI
interaction_avg_model_coef[1:2,survey_unit := "AI"][1:2,coef_true := coef][1:2,predictor := c("Relative fishing catch","Minimum temperature")]

#Other summed_tonnes_scaled_byreg
interaction_avg_model_coef[3:33,survey_unit := substr(coef_name, 39, str_length(coef_name))][3:33,coef_true := interaction_avg_model_coef[1,coef_true]+coef][3:33,predictor := "Relative fishing catch"]

#Other yearly_min_bypoint_avg.scaledacrossall
interaction_avg_model_coef[34:64,survey_unit := substr(coef_name, str_length("survey_unit")+1, str_length(coef_name)-str_length("yearly_min_bypoint_avg.scaledacrossall")-1)][34:64,coef_true := interaction_avg_model_coef[2,coef_true]+coef][34:64,predictor := "Minimum temperature"]

#reorder temperature and fishing
interaction_avg_model_coef[,predictor := factor(predictor, levels = c("Minimum temperature","Relative fishing catch"))]

#link for full survey name
interaction_avg_model_coef <- color_link[interaction_avg_model_coef, on = "survey_unit"]

#reorder by survey_unit
interaction_avg_model_coef[,Survey_Name_Season := reorder(Survey_Name_Season, survey_unit, decreasing = TRUE)]

Warning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NAWarning: argument is not numeric or logical: returning NA

#mark significance
interaction_avg_model_coef[,Significant := ifelse((coef_true-se > 0 & coef_true+se > 0) | (coef_true-se < 0 & coef_true+se < 0),T,F)]


#Plot both
sbt_fishing_avg_model_coef <- ggplot() +
  geom_hline(yintercept = 0, linetype = "dotted") +
geom_point(data = interaction_avg_model_coef, aes(x = Survey_Name_Season, y = coef_true, color = Significant)) +  
geom_errorbar(data = interaction_avg_model_coef, aes(x = Survey_Name_Season, ymin = coef_true-se, ymax = coef_true+se, color = Significant), width = 0) +
  scale_color_manual(values = c("darkgrey","black")) +
facet_wrap(~predictor, scales = "free_x") +
scale_x_discrete(limits = rev) +
  labs(y = "Coefficient", x = "") +
coord_flip() +
theme_classic()
  

Plot all other coefficients in averaged model

model_avg_values.nonfishortemp <- model_avg_values[c(2,3,4,100:104),]

#mark significance
model_avg_values.nonfishortemp[,Significant := ifelse((coef-se > 0 & coef+se > 0) | (coef-se < 0 & coef+se < 0),T,F)]

#more helpful names for variables
model_avg_values.nonfishortemp[,Variable := c("Area","Julian day","Species count","Number of tows","Latitude range","Depth","Latitude","Depth range")]

#make factor with order
model_avg_values.nonfishortemp[,Variable := factor(Variable, levels = c("Area","Julian day","Species count","Number of tows","Depth","Depth range","Latitude","Latitude range"))]


#plot
all_avg_model_coef <- ggplot() + 
geom_point(data = model_avg_values.nonfishortemp, aes(x = Variable, y = coef, color = Significant)) +  
geom_errorbar(data = model_avg_values.nonfishortemp, aes(x = Variable, ymin = coef-se, ymax = coef+se, color = Significant), width = 0) +
  scale_color_manual(values = c("darkgrey","black")) +
geom_hline(yintercept = 0) +
scale_x_discrete(limits = rev) +
  labs(y = "Coefficient", x = "\n\n\n") +
coord_flip() +
theme_classic()

#plot
all_but_lat_avg_model_coef <- ggplot() + 
geom_point(data = model_avg_values.nonfishortemp[Variable == "Latitude"], aes(x = Variable, y = coef, color = Significant)) +  
geom_errorbar(data = model_avg_values.nonfishortemp[Variable == "Latitude"], aes(x = Variable, ymin = coef-se, ymax = coef+se, color = Significant), width = 0) +
  scale_color_manual(values = c("darkgrey","black")) +
geom_hline(yintercept = 0) +
scale_x_discrete(limits = rev) +
  labs(y = "Coefficient", x = "Model variable") +
coord_flip() +
theme_classic()

#plot
lat_avg_model_coef <- ggplot() + 
geom_point(data = model_avg_values.nonfishortemp[Variable != "Latitude"], aes(x = Variable, y = coef, color = Significant)) +  
geom_errorbar(data = model_avg_values.nonfishortemp[Variable != "Latitude"], aes(x = Variable, ymin = coef-se, ymax = coef+se, color = Significant), width = 0) +
  scale_color_manual(values = c("darkgrey","black")) +
geom_hline(yintercept = 0) +
scale_x_discrete(limits = rev) +
  labs(y = "Coefficient", x = "Model variable") +
coord_flip() +
theme_classic()

#merge into single plot

model_coef_summary_sbt_jaccard <- cowplot::plot_grid(sbt_fishing_avg_model_coef+theme(legend.position = "null", axis.title.x = element_blank()),all_avg_model_coef+theme(legend.position = "null"), ncol = 1, labels = c("  a.                                                      b.","        c."), label_y = 0.99, rel_heights = c(3,1))

ggsave(model_coef_summary_sbt_jaccard, path = here::here("figures"),filename = "model_coef_summary_sbt_jaccard.jpg", height = 6.5, width = 8, unit = "in")